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Abstract. 

We investigate the energy transport in a one-dimensional lattice of oscillators with 
a harmonic nearest neighbor coupling and a harmonic plus quartic on-site potential. 
As numerically observed for particular coupling parameters before, and confirmed by 
our study, such chains satisfy Fourier's law: a chain of length N coupled to thermal 
reservoirs at both ends has an average steady state energy current proportional to 
1/N. On the theoretical level we employ the Peierls transport equation for phonons 
and note that beyond a mere exchange of labels it admits nondegenerate phonon 
collisions. These collisions are responsible for a finite heat conductivity. The predic- 
tions of kinetic theory are compared with molecular dynamics simulations. In the 
range of weak anharmonicity, respectively low temperatures, reasonable agreement 
is observed. 



1 Introduction 



In their seminal work of 1955, Fermi, Pasta, and Ulam ,1J investigate the relax- 
ation to equilibrium for a chain of coupled anharmonic oscillators, by exploiting the 
then newly available electronic computational devices. We refer to the informative 
memorial volume [2] . Related to their study is the issue of energy transport along 
the chain, which is modelled by coupling a chain of length N to thermal reservoirs 
at both ends. With increased computational power at hand such studies have been 
revived and carried through in considerable detail. Excellent reviews are in print 
[311] and here we only highlight, somewhat crudely, the main findings: 

(i) There are chains for which Fourier's law holds, in the sense that the steady state 
energy flux j e = 1/N. 

ii) There are chains with anomalous heat conduction, for which the energy flux ex- 
hibits a power law dependence on iV which differs from 1/N. 

iii) If the interaction depends only on the relative displacements, anomalous heat 
conduction seems to be the rule. 

iv) In some models the conductivity depends on the details of the coupling to the 
thermal reservoirs. 

While the amount of data available is impressive, it is generally agreed that 
there is very little theory which would serve as a guideline. The harmonic chain 
can be solved exactly with the result that j e is independent of N jjjj. For the 
harmonic chain with random masses the transport for a given wave number k is 
proportional to e~ 7 ^ k ' N , j(k) > 0, with 7 (A;) — > as k — * 0. Thus the average 
energy current depends on the precise spectral statistics of the thermal reservoir 
jHl U\ ■ For anharmonic chains there are attempts to predict the exponent for the 
anomalous heat conduction through mode-coupling theory E] . 

In our note we follow the strategy of Peierls, who argues that in case of weak 
nonlinearity one can use a Boltzmann type transport equation for the computation 
of the thermal conductivity. For anharmonic crystals in three dimensions phonon 
kinetic theory is well supported through theory ^0J El El arid also experimen- 
tally fOlEI]. Whether kinetic theory is applicable to a weakly anharmonic chain is 
somewhat tentative. There is no difficulty in writing down the appropriate trans- 
port equation. Its collision term ensures energy and momentum conservation in a 
phonon collision. To progress further an analysis of the solution manifold to both 
conservation laws becomes necessary. Firstly one has the trivial solutions in which 
the two colliding phonons merely exchange their label. For the label exchanging 
solutions the collision operator vanishes. However, as we first learned from Lefevere 
and Schenkel there is in addition a non-perturbative solution, which leads to 
nondegenerate collisions. As will be explained in more detail, with this input kinetic 
theory predicts a finite, non-zero, thermal conductivity even for a chain. 

As we learned later on, Pereverzev ^3] has already applied kinetic theory to 
the FPU /3-lattice, for which the potential energy depends only on the relative 
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displacements. For this model he obtains the non-perturbative solution and argues, 
based on the linearized transport equation, that the energy current correlation has 
a power law decay as t -3 / 5 . Thus the heat transport is anomalous with j e = N 2 ^ 5 . 
In our contribution we will be concerned with regular transport only. 

To keep the model as simple as possible we will study the harmonic chain with a 
small quartic on-site potential. Equivalently, one may consider a fixed anharmonicity 
but "low" temperatures. We regard the thermal conductivity to be defined through 
the Green-Kubo formula and work out the predictions of kinetic theory in fair detail. 
To explicitly compute the conductivity one has to invert the linearized collision 
operator. This is not completely straightforward and we have to be satisfied with 
some estimates, which however turn out to be sufficient for our purposes. The 
predictions of kinetic theory will be compared with molecular dynamics simulations. 
In fact, kinetic theory does rather well, with a range of validity larger than expected 
on the basis of purely theoretical arguments. 



2 The Green-Kubo formula, scaling properties 

The anharmonic chain is governed by the Hamiltonian 

H = Efe 2 + 5 w o?? - Sufaq+i + \Xqj} ■ (2.1) 

jez 

Here qj is the deviation from the rest position and pj the momentum of the j-th 
particle. We choose units such that the mass of a particle equals 1. u and 5u 
characterize the harmonic on-site and nearest neighbor interaction, respectively. 
ujq > and it has the dimensions of a frequency. To have a stable harmonic part of 
H we require 

0<<5< i. (2.2) 

In the border case 8 = ^, the harmonic part can be written as Ysji^o/tyiOj+i ~ Qj) 2 
and thus depends only on the relative displacement. A > is the strength of the 
quartic on-site potential. 

The particular case A = 1, 8 = ~ is studied in great detail in jTUJ , in which case 
kinetic theory is applicable at low temperatures. For our purposes it is of importance 
to add the extra parameter 5, since it is retained in the kinetic limit. Thereby one 
can compare the theoretical predictions with molecular dynamics simulations in 
their ^-dependence. 

To (j2.1j) we associate the local energy 

H 3 = \v) + + \Mj ~ l^ofe-i* + qjq j+ x) ■ (2.3) 

Then 

—Hj = —Jjj+i + Jj-ij (2.4) 
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with the energy current across the bond j to j + 1 given by 

J i,i+i = -fa 2 { Pj q j+1 - Pj+iQj) ■ ( 2 -5) 

The equilibrium state for ()2.1|) is Z~ x exp[— H/T] with T denoting temperature. 
Equilibrium expectations are denoted by (■)^ 1 ^ T \ We define the total energy current 
correlation through 

C(t;T,u Q ,5,X) = ^(4,+i(0)J ,iW) (1/T) - (2-6) 

By stationarity of (-)^, C(t) = C(—t). According to Green-Kubo the thermal 
conductivity at temperature T is defined by 

POO 

K (T) = T~ 2 / dtC(t-T). (2.7) 
Jo 

Regular transport in the sense of Fourier's law requires that < k(T) < oo. 

C(t) does not depend on all of its parameters separately. To find the dependence 
out we transform to new variables as 

<h(t) = l<lj( at ) , Pj(t) = aiPj(at) . (2.8) 
Then q(t), p(t) are solutions to Hamilton's equations of motion for the Hamiltonian 

H(q,p) = J2^P 2 J + ^o<?f - a 2 5uj 2 q 3 q J+1 + a 2 ^ 2 ^} . (2.9) 

In addition, 

H(q 1 p) = (a 1 y 2 H(q,p). (2.10) 

Therefore, 

C(t; T, u , 5, A) = (a-f 2 )- 2 a- A C(t/a; a 2 ^ 2 T, au , 5, a 2 7 " 2 A) , (2.11) 

which by (|2.7J) implies for the conductivity 

k(T, uj , 5, A) = a _1 fi:(a 2 7 2 T, aco , S, a 2 7~ 2 A) . (2-12) 

Setting auo = 1, a 2 7 ~ 2 A = 1 yields the scaling form 

k(T, luq, 5, A) = uqE,(uq XT, 5) . (2-13) 

In our molecular dynamics simulations, see Sect. 03 we set ujq = l/y/5 and A = 1. 
On the other hand, for the kinetic limit the natural choice is uq = 1 and fixed inverse 
temperature (3. Inserting a = \f5, 7 = (£/3T) _1//2 in 1)2.12)1 one arrives at 

T 2 k(T, * § , i) = IJ-^PTfKip- 1 , 1, 6, 6 2 f3T) . (2.14) 

Thus the limit T — > on the left corresponds to the limit A —>■ of A 2 /t(/5 _1 , 1, 5, A) 
on the right hand side, which will be studied in the following section. 
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3 Energy current-current correlation in the ki- 
netic limit 

In this section we set ujq — 1. To study the kinetic limit it is convenient to switch 
to Fourier space. Let T = [— |, \\ be the first Brillouin zone of the lattice dual to Z. 
For / : Z -> R we set 

f(k) = J2e- l2nkl f), keT, (3.1) 



with the inverse 



We decompose 



fi = f dke l2 ^f(k) . (3.2) 



H = H hax + ±\J2<lj- (3-3) 
The harmonic part iJhar has the dispersion relation 

uj{k) = {l-25 cos(2ttA;)) 1/2 . (3.4) 
PjiQji 3 e ^5 are concatenated into a single complex- valued field through 

a(fc) = ±={^(k)q{k) + i-±—ftk)) . (3.5) 

Notationally it will be convenient to also define 

a(k)* = a{k, 1) , a(Jfe) = a(k, -1) . (3.6) 
Then the equations of motion read 



d_ 

dt 



a(k,t) — —iuj(k)a(k,t) — iX 2_. \ 



fl,0'2iO'3=±l 



dkidk 2 dk 3 

T 3 
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(16cj(A;)cj(A;i)cj(A; 2 )cj(A;3))" 1/2 5(A; + <7ifci + ^2 + ^3) JJ a(%, o-j, i) . (3.7) 

3=1 

In these new variables the harmonic part of the energy becomes 

H har = ! dku(k)a(k)*a(k) (3.8) 
J 

and the total current, J tot = £V Jjj+i, becomes 

J tot = S f dksm(2nk)a(k)*a(k) = [ dk(^-Lu(k))Lu(k)a(k)*a(k) . (3.9) 



2n l T dk 



We plan to study the energy current correlation (|2.6J) in the limit of small A. 
Since in equilibrium (Jo,i(^))^ = 0> ^ holds 

E^iW-W^ = lim -(J ,iW) (ftr) , (3.10) 

where (•)^ 3 > r ) refers to expectation with respect to the perturbed equilibrium measure 
Z^ 1 exp[— (3H + t Jtot]> ^tot the total current, in the infinite volume limit. For small 
A we can ignore the quartic on-site potential in the average (■)^ l3 ' T \ Thus (•}^' r ) is 
replaced by the Gaussian measure (-)^ rm , which is uniquely characterized through 
its covariance 

(a(*)*a(*0>2£m = ^ k ~ *O0M*0 - r5sin(2vrA;))- 1 , 

<«(*)>£L = °> (fl(*M* / )>£L = o- (3- 11 ) 

The state (■)^ m is not invariant under the mechanical time evolution and we have to 
understand how for small A such a non equilibrium measure evolves in time. Instead 
of (-)farm ^ us consider as initial state an arbitrary translation invariant Gaussian 
measure with covariance 

(a(AO) = 0, (a(AOa(fcO)o = , 

(a(k)*a(V)) Q = 6(k - k')W(k) , (3.12) 
compare with (|3.11|) . Its two-point function at time t is given by 

(a(k, t)*a{k', t)) = (a(k)*a(k')) t = 5{k - k')W x {k, t) , (3.13) 
which defines W x (k,t). In particular, by (J3.9)) . 

(J 0jl (t)} = <5 / dksm(2nk)W\k,t) . (3.14) 



For A = 0, one has W °(k,t) = W(k). For small A the time variation of W x (k,t) 
is on the time-scale A~ 2 . Thus one expects that the following limit exists, 

\imW x (k,\-H) = W(k,t) , (3.15) 

A — >0 

and that the limit Wigner function W(k, t) evolves according to the spatially homo- 
geneous Boltzmann equation 



d_ 

dt 



W(k,t) = 12n / dkidk 2 dk 3 {lQujuJiUj 2 u 3 ) 1 



""1,0-2,(73 = 

5(lu + o\LO\ + a 2 uo 2 + a 3 uj 3 )5(k + o\ki + a 2 k 2 + a 3 k 3 ) 
[WiW 2 W 3 + W{a x W 2 W 3 + a 2 W l W 3 + o 3 W x W 2 )\ (3.16) 
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with initial conditions W(k,0) = W(k). Here, and later on, we use the shorthand 
u = u(k), ojj = uj(kj), W = W(k), Wj = W(kj), j = 1,2,3. In Appendix A we 
explain the second order diagrammatic expansion in A, mostly to make sure that 
the collision strength is correct. 

By the argument in Appendix 18.1 of energy and momentum conservation 
in ()3.16|) can be satisfied only if 

3 

i.e., only for phonon number conserving collisions. Hence ()3.16|) simplifies to 
d 9n f 

-W(k, t) = — dkidk2dk^,{ujuJiuj2UJ-i)~ 1 5{uj + lo\ — uj 2 — 0Jz)8(k + k\ — k 2 — k 3 ) 



dt y ' ' 4 . 

[WiW 2 ^ 3 + W{-W X W 2 - W x Wz + W 2 W 3 )] 
= C(W{t)){k). (3.18) 

Clearly, the equilibrium Wigner function 

Wp(k) = -4tt (3-19) 

is a stationary solution for ()3.18|) . Since according to ()3.10|) and (j3.11|) only small 
deviations from equilibrium are needed, we linearize the collision operator as 

C{W P + (W p ) 2 f) = -{3~ A Lf + O(f) . (3.20) 

Then 

97T f 

Lf(k) — — I dkidk 2 dk 3 (uJUJiUJ2UJ 3 )^ 2 5(to + U)\ — cu 2 — UJ3) 

5(k + h-k 2 - k 3 ) (f(k) + f(h) - f(k 2 ) - f(k 3 )) . (3.21) 

For this particular choice of linearization L = L* in L 2 (T, dk). 
Let Vf(k) = Wp(k)f(k). Then 

C(W P + f) = -Af + O(f) = -(3- 4 LV~ 2 f + 0{f) . (3.22) 

Let g{k) = sin(27rfc) and let (•, •) be the scalar product in L 2 (T, dk). Then, with 
W T (k,t) denoting the solution to ()3.18|) for the initial condition W(k) = \Bu{k) — 

tS sin(27rA;)) , one has 

lim C{\~ 2 t- p~\ 1, 6, X) = lim -5{g, WJt)) 

A^O r^O T 

= 5 2 (g,e- M V 2 g) 

= 5 2 (g,Ve-^ y - lLV ' lt Vg). (3.23) 
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Integrating over time one concludes that 

lim X 2 P' 2 k{(3-\ 1, 5, A) = 5 2 f3 i (g, V^^LV'^Vg) 

= S 2 ^(g, V 2 L^V 2 g) 

= 6 2 (tu- 2 g,L- 1 tu- 2 g). (3.24) 
Combining with (|2.14j) yields as the final result 

KmT\(r, -J=, 5, 1) = 5- 5/2 (uj- 2 g, L~ l uj~ 2 g) . (3.25) 

The inner product on the right depends only on 5. The prefactor is slightly mis- 
leading, since L is proportional to Hence 

8- s ' 2 (u- 2 g, L- x u- 2 g) = cT 3 / 2 c(5) (3.26) 

with < c(5) < oo. 

()3.24|) can rephrased as the asymptotics of the scaling function from (j2.13jh 

lim x 2 E(x, 5) = 5 2 {u~ 2 g, L^u^g) . (3.27) 

In other words 

k(T, lo , 5, A) S (cu ) 9 (\T)- 2 5 2 (cu- 2 g, L'^g) (3.28) 

for small cUq XT. Our argument provides no indication over which range ()3.28|) 
is a valid approximation. In fact, even the claim ()3.25|) . ()3.26|) is tentative. The 
diagrammatic expansion from Appendix A relies on the separation into leading and 
subleading diagrams. For a three-dimensional lattice such a separation is convincing 
and can be checked for special diagrams. The oscillatory time integrals for the chain 
have a slower decay and the rough estimates used so far are not sufficient to justify 
the separation into leading and subleading diagrams, which we have assumed here. 
On the other hand, there could very well be cancellations which are difficult to 
access through the expansion in Feynman diagrams. Also we need the validity of 
the kinetic equation only close to thermal equilibrium. In view of this situation a 
molecular dynamics simulation is in demand. 

4 Thermal conductivity from the one-dimensional 
Boltzmann equation 

The linearized collision operator L of (|3.21|) corresponds to the quadratic form 



1 9n f 

(f, Lf) = -:-7-\ dk 1 dk 2 dk 3 dk 4 ,(uj 1 uj2UJ 3 U4)~ 2 S(uj 1 + io 2 - to 3 - u A ) 
4 4 J T 4 

S(h + h-k 3 - h)(h + f 2 -h- Uf ■ (4.1) 
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Figure 1: The non-perturbative solution h(ki) k 3 , S) as a function of k\ for fixed 
k 3 and 5. On the left, 5 = 0.4, and on the right, 5 = 0.5. For each value three 
solutions are plotted, corresponding to k 3 = 0.02 (black), 0.25 (dark grey), and 0.5 
(light grey). We plot T x T in the extended zone scheme, the dashed lines are the 
boundaries of a unit cell. For 5 < 0.5 the non-perturbative solution is smooth, while 
for S = 0.5 there are cusp singularities. 



We use momentum conservation to integrate over fc 4 . For energy conservation we 
thus need the solutions to 

w(fci) + u{k 2 ) - uj{k 3 ) - u(k! + k 2 - k 3 ) = . (4.2) 

The obvious solutions are 

ki = k 3 and k 2 = k 3 . (4.3) 

Expanding (J4.2)) to first order in 6, as noticed in ^3], there is yet a further, 
non-perturbative solution given by 



ki + k 2 = - modulo 1 for all k 3 . 



This suggests to write, in general, 



k 2 = h(ki; k 3 ) . 



(4.4) 



(4.5) 



Indeed, for every k 3 G T there exists a unique function h(-; k 3 ) : T — > T, which is 
continuous, one-to-one, and satisfies 



uj(ki) + u(h(ki, k 3 )) - u(k 3 ) - w(ki + h(fa; k 3 ) - k 3 ) = . 



(4.6) 
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h is called the non-perturbative solution. In Fig.^we display a few non-perturbative 
solutions at 5 = 0.4 and 5 = 0.5 for three values of k 3 . For small 5 one finds 

h(h; k 3 ) = ]--k 1 -8^-{ sm^TTh) + sin^TrA*)) + 0(8 2 ) , (4.7) 
A 271 

which reasonably well approximates the left hand of Fig. ^ 

Inserting the solutions (|4.3jl and ()4.5|) to energy conservation into (j4.1|) splits the 
linearized collision operator as the sum 

L = L cx + L npert . (4-8) 

By symmetry, for the label exchanging solution (/, L ex f) = for all /. Hence 
L ex = and energy conservation in (|4.1|) will be evaluated always at k 2 = h{k\] k 3 ). 

In phonon kinetic theory it is customary to distinguish between normal and 
umklapp processes. We choose the convention that kj G [— ~, ~], j = 1,...,4. 
Then a process is called normal if k\ + k 2 — k 3 — k 4 = 0, while it is umklapp if 
fci + &2 — &3 — &4 = ±1. By this definition, the curves in Fig. ^are divided into a 
normal piece and an umklapp piece. For example, at 5 = 0, one has k\ + k 2 = ±1/2 
and k 3 + fc 4 = ±1/2. If the two terms have opposite sign, the collision process is 
normal and otherwise it is umklapp. In our context such a division looks artificial. 
In fact we will find that both, normal and umklapp, processes contribute to the 
thermal conductivity. 

We turn to the zero subspace of L, i.e., to solutions of Lf = in L 2 (T). By 
(|4.ip . clearly they must be collisional invariants in the sense that 

f(h) + f(h(k i; k 3 )) - f(k 3 ) - f(h + h(h; ks) -h) = (4.9) 

for all (hi, k 3 ) 6 T 2 . The obvious solutions are 

/(*) = !, f(k)=u(k). (4.10) 

We expect that there are no further solutions, but no proof is available, at present. 
This is an important issue, since speaking in general, the number of collisional 
invariants is the crucial information on the long-time behavior of a kinetic equation. 
At S = 0, h does not depend on k s and as a consequence the zero subspace of L 
becomes infinite-dimensional consisting of all /'s satisfying f(kx) + /(- — ki) = 0. 

We integrate in (|4.ip over fc 4 and k 2 . For the volume element with respect to k 2 
we need 

A (u(h) + u(k 2 ) - u(k 3 ) - uj(h + k 2 - k 3 )) \ k2=h{kl . M) (4.11) 
= (2n5) (lu^)- 1 sm(2nk 2 ) - uo(k x + k 2 - A; 3 )~ 1 sin(27r(A; 1 ± k 2 - k 3 ))) 
Hence 

- ^— l / jl jl /, ,/u \, j v. \. j v. \, ,n„ i l l ^^ _2 



k2=h{k\\kz) 



(/, Lf) = -rr:{^5)- 1 / eZMfcaMfciM^MfcaMfci + h ~ h)) 

ID Jj2 

Itoih)- 1 sm(2nk 2 ) - u(k x + k 2 - k 3 y l sin(27r(A;i + k 2 - h))]' 1 
(f(h) + f(k 2 ) - f(k 3 ) - f(h + k 2 - k 3 ))\ =h{ki . k3) (4.12) 
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and L carries an explicit pref actor 8 , as claimed in 1)3.26)1 . 

To obtain the thermal conductivity (in the kinetic regime) one has to invert L, 
which can be achieved only numerically and which is not completely straightforward 
because of the constraint due to energy conservation. However for small 5, say up 
to 5 = 0.35, more modest means already suffice. By Jensen's inequality one has 



{f,L- 1 f)>{fJ) 2 /{f,Lf) (4.13) 
with f(k) = u(k)~ 2 sin(27rfc). For 5 = one obtains 

-1/2 , 
' ^ .,,■>,,< 2 I 

"1/2 



(fj) 2 = ( [ 1/2 (sm2nk) 2 dk) 2 = ] (4.14) 

</-l/2 4 



and 

-1/2 rl/2 



(f, Lf) = r 1 J- / dh [ dk 3 \ sin(27rA;i) - sin(27r£;3 

61 </-l/2 J -1/2 

, 2 



(2sin(27rfci) - 2sin(27rfc 3 ))' 



= 5- l - 2 . (4.15) 

7T 

Combing ()3~23)) and (gHH) yields 

lim T 2 k(T, -^=, 5, 1) > 5- 3/2 — 5" 3/2 0.27 (4.16) 

for small 5. 

For numerical inversion of L at 5 = we expand in a basis of the form sin((2n + 
l)2nk), n = 0, 1, . . ., and, instead of c(0) = 0.27, obtain the prefactor in (j3.26j) as 
c(0) = 0.275637 with a stable value over the range n = 5, . . . , 32. As can be seen 
from Fig. ^ even at 5 = 0.4 the approximate solution ()4.7|) is rather accurate. 
Therefore c(S) is expected to depend only weakly on 5. Our claim is supported by 
the lower bound 

c(5) = 5~ 1 (u~ 2 g, L~ 1 u~ 2 g) > (u^g, ou^g) 2 /S^g, Lu' l g) , (4.17) 

g(k) = sin(27rfc). For this particular choice of the variational function the singular 
denominator in (|4.12|) is cancelled exactly and the numerical integration, using the 
true non-perturbative solution, becomes routine. The lower bound to c(5) drops 
from 0.27 at 5 = to 0.2 at 5 = 0.35. More details on the linearized collision 
operator can be found in [T7j . 

5 Thermal conductivity from molecular dynamics 

Kinetic theory is expected to be valid for small dimensionless coupling Uq 4 \T and 
large system size. How small a coupling and how large a system size can be explored 
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Figure 2: The relation between the average current J/N and the temperature gra- 
dient, VT, for T = 0.4,5 = 0.2, N = 200. The straight line is Fourier's law as 
obtained by a fit with gradient k. 

only through molecular dynamics simulations. To this end, we compute thermal 
conductivities for various parameters of the anharmonic chain (j2.1j) . We set A = 1, 
ujq — 1/ VS. Then the infinite volume conductivity depends only on 5, T and kinetic 
theory becomes valid in the limit T — > 0, compare with Sects. El and El In particular, 

T 2 k(T) ~ <T 3/2 c(5) (5.1) 

for small T with c(0) = 0.28 and c(5) slowly dropping to smaller values as 5 is 
increased. 

To numerically determine the thermal conductivity we take a chain of finite 
length, N, and attach thermostats at both ends. We use free boundary conditions, 
which means = qj at the boundaries, but the physical results are insensitive to 
this particular choice of boundary conditions. We adopt deterministic thermostats 
which generalize those of Nose-Hoover [IBllISj, arid follow the methods used for the 
chain when 5=1/2 fHlEHj- The non-equilibrium steady state is achieved by inte- 
grating the equations of motion and physical observables are measured by averaging 
over time after waiting for a sufficiently long equilibration time span. The inte- 
gration is performed numerically using standard algorithms, such as Runge-Kutta, 
with time steps of 0.001 ~ 0.04. 10 7 ~ 10 10 samples are taken to obtain the average 
values of physical observables. The results do not depend on the time step size. A 
non-equilibrium steady state has as its parameters the boundary temperatures, N, 
and 5. The local temperature is defined through the relation Tj = (p 2 )- 
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Figure 3: The autocorrelation function (Jo,i(0) Jj t j + i(t)) at t = 0,20,40,60 for T = 
OA, 5 = 1/3. The peaks move away at a constant rate, from which the speed of 
phonons is extracted. 



In principle, a single non-equilibrium state suffices to obtain k. To determine k 
with enough accuracy we employ a more elaborate procedure. We choose boundary 
conditions with varying temperature gradients, VT, but with the same temperature 
close to the midpoint. We can then check that the average current is linear in VT 
and obtain k for given T, S, and N, as illustrated in Fig. El The temperature gradient 
needs to be computed away from the boundaries because of jumps in the temperature 
at both ends, see Sect. |B] for a discussion. Since local energy is conserved and there 
are no internal heat sinks or sources, the average current is constant throughout the 
system. We increase N and examine if the bulk limit is reached to finally extract k 
for given T, 5. Bulk behavior had already been observed when 5=1/2 [2011^], and 
we do so also in the present study. 

Naively, it would seem that weak coupling physics should be easy to simulate. 
The computational difficulties arise because the mean free path inevitably becomes 
large and the relaxation towards the steady state becomes slow, which demands to 
carry out the simulation for large system sizes and sufficiently long times. These 
factors limit the range of accessible parameters values. To elucidate these points 
and also to understand the kinetic theory aspects of this system better, we analyze 
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Figure 4: Average phonon velocity, v, for T = 0.1 (□), 0.4 (O), 4(A), 10(x). The 
solid line represents the kinetic approximation (J5.4|) . 



the statistical properties of the system in more detail. To judge the required system 
size we have to estimate the mean free path A. It is not a sharply defined quantity. 
Following Ziman ^U] we set 

A = k/(C v v), (5.2) 

where C v is the specific heat and v the average speed of phonons. In the harmonic 
approximation one obtains C v = 1. Kinetically the equilibrium phonon number 
density equals T/u(k), which suggests to set 



v = ( J (^-fTuj^dk I J Tu^dk 

with the following expansion for small 5, 

^ A 9 



1/2 

(5.3) 



2 [} + f/ + ° {64) ) ■ (5 ' 4) 

The speed of phonon propagation can also be measured directly from the time 
and space dependences of the autocorrelation function (Jo,i(0) Jj t j + i(t)) in thermal 
equilibrium 1X0] - The velocity of the peaks in the correlation function is equated 
with the average phonon velocity relevant for thermal transport. An example is 
shown in Fig. El The measured velocities can be compared to the kinetic result 
in (j5.4|) . which is done in Fig. |3J Perfect agreement is not expected for a number 
of reasons. First, there is no unique definition of the average speed of phonons so 
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Figure 5: T 2 k(T, l/Vd, 5, 1) compared against the kinetic theory computation, 
0.28<5~ 3/2 (straight line). The data points are for T = 0.1 (□), 0.4 (O), 4(A). 

that there is no guarantee that the average ()5.3|) is precisely what we measure in the 
molecular dynamics simulations as in Fig. El Furthermore, the kinetic approximation 
for v does not include the effects of anharmonicities and should strictly hold only in 
the limit T — > 0. Given these constraints, the agreement between the simple formula 
()5.3j) for v and the simulation results in Fig. 0] is fairly satisfactory. 
From the above discussion and using ()5.1|) we obtain 



for small 5, T. From the measured values of k, v (see Figs. |U and EJ), A can be 
determined according to (j5.2j) . A ranges then from A ~ 1 at T = 4, 5 = 0.3 to 
A ~ 4000 at T = 0.1, 5 = 0.08. In the simulations the system size is varied up to a 
few thousand depending on the parameters. Therefore, in our simulations, we have 
been able to achieve N > A in the parameter range probed here. When A and N are 
of the same order, it is not clear a priori if the bulk limit has been reached. In the 
subsequent section, we provide an argument that we might still be able to estimate 
the conductivity in such cases, even when N < A. 

One distinctive qualitative feature of kinetic theory is the leading <5 -3 / 2 depen- 
dence of T 2 k(T), see (|5.1|) . To compare the molecular dynamics simulation results 
with kinetic theory, we need to keep in mind that the agreement should hold only 
when T is small and the power law becomes exact when 5 —>■ 0. In Fig. the 
comparison is made and we note a surprisingly good agreement, which improves as 
the temperature is lowered, as to be expected. 



A ^ 0.38 cT 2 T~ 2 



(5.5) 
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6 Steady state temperature profile 

A second quantity of physical interest is the average temperature profile. If dif- 
ference in the boundary temperatures, AT = T_ — T + , is small, then the profile 
is approximately linear. For larger AT the temperature dependence of k will be 
seen. As AT is further increased, local equilibrium will break down eventually [22J. 
We begin with working out the temperature profile as predicted by the transport 
equation. 

The space-time dependent version of the Boltzmann equation (j3.18|) reads 

-W(x, k, t) + —u/(k)—W(x, k, t) = C(W(X, ; t))(k) , (6.1) 

Ot Z7T OX 

cu'(k) = -^to(k), where the collision operator is local in (x, t), i.e., it acts only on the 
wave number k. Under (J6.1|) the phonon number density 



T 



n(x,t) = I W(x,k,t)dk (6.2) 

and the energy density 

e(x,t)= I cu(k)W(x,k,t)dk (6.3) 



T 



are locally conserved. The former one we regard as spurious, since it has no analogue 
on the microscopic level. 

Following the standard hydrodynamic scheme, since there are no convective 
terms, the long time behavior of (|6.1j) is thus dictated by the solution of the coupled 
nonlinear diffusion equations 

9M d -D(a, ffliM. (6.4) 



dt \e J dx dx \P 

Here a, (3 are "chemical potentials" labelling the stationary solutions, W a ^(k), of 
(ffi~Tj) as 

W a , p (k) = (puj(k) + a)' 1 (6.5) 

for (a, (3) G V, where V = {a, P\u(0)P > -a for (3 > and u(\)f3 > -a for (3 < 0}. 
Then V 3 (a, (3) i— > (n, e) with 



n(a, p) = W a> p(k)dk , e(a, p) = u{k)W a>p {k)dk . (6.6) 

In (|6.4j) we insert the inverse function as defined on (M+) 2 , which is uniquely specified 
because of convexity. Secondly, D(a,P) is the 2x2 matrix of Onsager coefficients 
as given through a Green-Kubo formula analogous to (|3.24|) . Following in spirit the 
arguments from ^2] one obtains 

D(a, P) = (2n)- 2 (co'(W a ^ 2 Q , I (1 u)u'{W a ^f) (6.7) 
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Figure 6: Steady state correlations (piPj) — TjSij for boundary temperatures T_ = 
0.8, T + = 1.2, 8 = 1/3, and N = 40. 

with the linearized collision operator 

Lf(k) — — I dk 1 dk2dk 3 (touJiU!2t03)~ 1 5(cu + — u; 2 — cj 3 ) 

4 Jj3 

5(k + h-k 2 - k^W^W^ik^W^^W^ih) 
{f(k)+f(k 1 )-f(k 2 )-f(k 3 )). (6.8) 

At a = 0, D(0,/3) = D ee is independent of (3, which results in an important 
simplification. Let us consider the steady state problem for ()6.4|) in the slab [0, 1] 
with boundary conditions a(0) = 0, a(l) = 0, /5(0) = (3{1) = (3 + . Then the 
solution is given by 

a(x) = 0, (3{x) = (3-{l -x) +(3+x, < x < 1 . (6.9) 

In particular the steady state energy flux is 

3e = D ee (P + -/3-). (6.10) 

For a small temperature difference, /3_ = T _1 , j3 + = (T + AT)~ X , one has in 
approximation 

Je = -D ee T- 2 AT, (6.11) 

where D ee = (u~ 2 g , L^u^g) in agreement with (|3.25|) . 

In molecular dynamics simulations the two ends of the chain, j = 1 and j = N, 
are coupled to thermal reservoirs with AT/T of order 0.4 or less. The local value 
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Figure 7: The temperature profile for 5 = 0.2, boundary thermostat temperatures 
T_ = 0.8, T + = 1.2, and N = 400. The dashed line is a fit with a linear l/3}-profile. 

of a is extracted from the simulation. A simple test is provided by the momentum 
covariance. If locally the Wigner function has the form W a ^(k), then 

(PiPj) = / (/MA;) + a)~ l u{k) cos(27r(i - j)k)dk , (6.12) 

which reduces to (piPj) = [3~ l 5ij for a = 0. Expanding ()6.12|) in a, 5 results in 
(pjPj+i) j (p?) ~ —a5/2. Numerically, the steady state momentum correlation (piPj) 
is indeed strongly peaked at i = j and decays rapidly. The simulations therefore 
indicate that a is small compared to u and hence one can safely set a = 0. An 
example for off-diagonal correlations obtained from molecular dynamics simulations 
is shown in Fig. Inland, in this case, \(pjPj + i) \ f(pf) ~ 0.02. This quantity is close to 
but not quite zero. This may be due to a being small, but not quite zero, or due to 
finite size corrections to local equilibrium. 

In Fig. 0we display a numerically generated steady state profile. One notes that 
the profile lies slightly below the linear interpolation. Indeed, (16 .9|) claims that 1/3} 
is linear. The dashed line is the corresponding fit. The excellent agreement is a 
further confirmation for a being small. 

A generic temperature profile consists of three pieces: there are two boundary 
jumps of equal size and concentrated over a few lattice spacings and there is an, in 
essence, linear bulk piece. Only if the effective temperature difference is large, one 
observes deviations from linearity, as for example in Fig. [7| According to [H] , and 
as easily extended to the case under study, in the harmonic limit boundary jumps 
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would be exponentially localized and the bulk piece would be flat. The observed 
temperature profile deviates significantly from this harmonic limit for all data points 
from Fig. |5] in always having a nonzero slope. 

The equal size boundary jumps are easily understood in the context of Langevin 
reservoirs. The leftmost particle, j — 1, is then governed by the Langevin equation 

Qi = Pi , 

pi = -ul qi + 5u 2 q 2 - Xqf - 7 pa + {2 1 TJ) 1 ' 2 b{t) , (6.13) 

where 7 > is the friction coefficient and b(t) white noise such that (b(t)b(t')) = 
S(t — t'). A corresponding equation holds for the rightmost particle, j = N, with T_ 
replaced by T + . In the steady state (Hi) = 0, which implies 

7 (T_-Ti) = J(N), (6.14) 

where J(N) is the steady state current for chain length N. Similarly 

j(T N -T + ) = J(N). (6.15) 

In particular, the boundary jumps are equal. 

To understand the full structure of the steady state one may resort to (|6.1J1 
restricted to the slab [0, £]. Then the steady state Wigner function, W(x, k), satisfies 
dW/dt = with the following boundary conditions at the two ends, x — 0, £, 

W(0, k) = ^uj'(k)T_/uj(k) for k > , 

W(£, k) = — -u)'{h)T + /uj{k) for k < , (6.16) 

where we assumed thermal sources together with complete absorption. To obtain the 
steady state profile on this basis would need considerable numerical effort. Therefore 
we turn to a very much simplified model, which however retains the gross features 
of the steady state. 

Energy is transported to the right with velocity +1 and to the left with velocity 
— 1. Locally the velocity may switch its orientation randomly with rate 77. Then 
the mean free path is I/77 and, as we will see, Fourier's law holds with conductivity 
k = 1/ (277) . In the steady state the average energy density satisfies 

A/ + (s)=77(/4*) -/+(*)), 

~f-(x)=r,{f+(x)-f.(x)). (6.17) 

The local energy is f+(x) + f-(x) = T(x), which we identify with the local temper- 
ature. The energy current is then f+(x) — f-(x). Energy is injected at x = I and 
absorbed at x — 0, i.e. 

/_(*) = !, / + (0)=0, (6.18) 
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where partial absorption could readily be included. Comparing to the case /+ (0) = 1 
and f-(£) = 1, one concludes that the imposed right boundary temperature is T + = 2 
and correspondingly the imposed left boundary temperature T_ = 0. The solution 
to (fHT7jl . dSHD reads 

f+(x) = (l + riiy^x, f4x) = (l + 7i£)- 1 (rix + l), (6.19) 

which yields the temperature profile 

T(x) = f + (x) + = (1 + 7 1 e)- 1 (27 1 x + 1) (6.20) 

and the steady state current J{£) for slab length £ as 

J(i) = f + (x)-f4x) = (l + r ] e)- 1 . (6.21) 

Thus at both ends the boundary jump equals and the effective temperature 

difference is AT = (1 + r)£)~ l 27]£. Hence 

J(£) = -kAT/£. (6.22) 

Thus, even if £ r]^ 1 , the correct bulk current is extracted through ()6.22|) . Of 
course, as for fixed £ the mean free path increases, so does the relaxation time and 
longer simulation times would be needed. 

To come back to the molecular dynamics computation, the simulation time is 
sufficiently long so that the steady state is reached. For most of the data presented, 
we obtain the conductivity from lattices with N ^> A and are confident that they 
should be reliable. At T = 0.1 data with small 5 are computed for lattices with 
N ~ A. Given the reasoning within the simplified model and the tendency of the 
data for larger T, we believe that a reasonable estimate for the conductivities has 
been obtained even for these cases. 

7 Conclusions 

The Boltzmann-Peierls equation retains the exact dispersion relation and the type 
of nonlinearity of the Hamiltonian model. For example, had we considered a cubic 
nonlinearity, then Eq. ()3.18|) would be quadratic in W. For a nonlinearity which 
depends on the nearest neighbor relative displacements there would appear the addi- 
tional factor | FJ^ =1 (1 — exp(i27r/c :) '))| 2 in the collision rate. It is remarkable that with 
this input the qualitative features of the "low" temperature thermal conductivity 
are so well recovered. 

It seems to us that the Boltzmann-Peierls equation has never been tested in 
comparable precision before. The peak time for the experimental investigation of 
phonon thermal conductivity was in the late 50'ies and early 60'ies. A quantita- 
tive comparison with the theory was hampered from two sides: (1) The dispersion 
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relation and the anharmonicities of a given dielectric crystal are not so readily avail- 
able. (2) One needs considerable numerical effort to reliably obtain the thermal 
conductivity from the transport equation. Thus mostly one had to be satisfied with 
qualitative predictions, as for example the 1/T-dependence of the thermal conduc- 
tivity in the presence of only three-phonon collisions [T0J. More recently, molecular 
dynamics simulations have become available, for example see |23I and references 
therein. Compared to these more material science oriented studies, we achieve a 
much larger system size, due to one instead of three spatial dimensions, and we test 
the simulation data directly against the transport equation without further approx- 
imations. 
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A Appendix 

We follow Sect. 11 of jT^j. The initial measure is translation invariant and uniquely 
characterized by the covariance of ()3.12|) . We want to establish that 

W x (k, \~H) = W{k) + tC(W)(k) + 0(f) (A.l) 

for small A. 

By (|3.7j) the vertex weight is given by 

<f)(k, h, k 2 , h) = (IQuju^ujz) 1 ' 2 . (A.2) 

We use the expansion through Feynman diagrams as obtained from the iteration 
of (11.7) in [12 j . Since the Hamiltonian has a quartic nonlinearity, each interaction 
leads to a branching into 3 lines, compare with Fig. 03 of [T2*j . 

To order A we simply have S(q — p)W(q). The order A term is purely imaginary, 
thus vanishes, because to every diagram there is its complex conjugate, denoted by 
c.c. Thus we are left with the order A 2 . It has 8 ways of branching. For a given 
branching there are 15 Gaussian pairings and 8 possible orientations of the internal 
lines, which in total amounts to 960 diagrams. They will be divided into subleading 
and leading. 

(i) There are 144 diagrams of the type 
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We set t = t 2 — t\. Their sum is then 



36 V 5(q-p) [ dfadk 3 (f)(q,fa,q,k 3 ) 2 W 1 W(-a 2 q)W 3 (e iT ^ 1+ ' T ^ + c.c.). (A.3) 



<72=±1 

(ii) There are 144 diagrams of the type 
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Their sum is 



36 V 5(q-p) [ dk 1 dk 3 (p(q,k 1 ,q,k 3 ) 2 W 1 W(q)W 3 a 2 (e iTUj{q)(1+a2) + c.c.) . (A.4) 



C72=±l 

(iii) There are 288 diagrams of the type 
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Their sum is 

36 V S(q - p) [ dfadk 3 (f>{q, fa, q, fafW^ {q)W 3 a 2 {j™^ + ^ + c.c.) = . 

(A.5) 

In (i) and (ii) the terms independent of r cancel each other. The remaining 
terms are proportional to cos(2o;(g)r) and thus of order A 2 after time-integration. 
We are left with 384 leading diagrams. 
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(iv) The gain term results from 96 diagrams of the type 
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Their sum is 



6 / S(q —p) dk 1 dk 2 dk 3 (f)(q, ki, k 2 , k 3 ) 2 5(q + 0\kx + a 2 k 2 + o- 3 k 3 

a"l,0-2,CT3 = ±l 



WiW 2 W 3 (e iT{ - uj{q)+aiUl+a2U2+a ' AUl3 ' ) + c.c. 
(v) The loss term results from 288 diagrams of the type 



(A.6) 
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Their sum is 



6 S(q-P) J dkidk 2 dk 3 (f)(q,k 1 ,k 2 ,k 3 ) 2 S(q + aih + a 2 k 2 + a 3 k 3 ) 

W{q){a 1 W 2 W 3 + a 2 W x W 3 + a 3 W 1 W 2 )(e lT{u{q)+(7lU)1+ ' T2U)2+ ' T ^ ) + c.c.) .(A. 7) 



In (iv) and (v) we use that 

~\~ 2 t r t2 



lim A 2 



dU 



o 



dtl ( e M*2-ti) + cc ) = 2vrt(5(cu) 



(A.* 



o 



when integrated against a smooth, rapidly decreasing test function. By adding (iv) 
and (v) one obtains the collision operator from Eq. ()3.16|) with the prefactor 12n. 
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